Improved models of stellar core collapse and still no explosions: 

What is missing? 
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Two-dimensional hydrodynamic simulations of stellar core-collapse with and without rotation are 
presented which for the first time were performed by solving the Boltzmann equation for the neutrino 
transport including a state-of-the-art description of neutrino interactions. Although convection 
develops below the neutrinosphere and in the neutrino-heated region behind the supernova shock, 
the models do not explode. This suggests missing physics, possibly with respect to the nuclear 
equation of state and weak interactions in the subnuclear regime. However, it might also indicate a 
fundamental problem of the neutrino-driven explosion mechanism. 
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Despite of more than three decades of theoretical 
research and numerical modeling, the processes which 
cause the explosion of massive stars are still not under- 
stood. Current observational data of supernovae (SNe) 
do not provide direct information. Although neutrinos 
(y) and gravitational waves could yield such insight, the 
v events detected in connection with SN 1987A were 
too sparse to constrain the SN mechanism. Progress in 
our understanding of the complex phenomena in collaps- 
ing stars and nascent ( "proto-" ) neutron stars (PNSs) is 
therefore mainly based on hydrodynamic simulations. 

Stars more massive than about 10 solar masses (M ) 
develop an iron core in a sequence of nuclear burning 
stages. This iron core becomes gravitationally unstable 
when it reaches a mass close to its Chandrasekhar limit 
and collapses to a neutron star. A hydrodynamic shock 
forms when nuclear density is reached and the matter be- 
comes incompressible. There is general agreement, sup- 
ported by detailed numerical models, that this shock is 
not able to promptly cause a SN explosion. Instead, it 
suffers from severe energy losses by the photodisintegra- 
tion of iron nuclei to free nucleons. It finally stalls after 
having reached densities low enough that electron neutri- 
nos (y e ) can rapidly escape in a luminous outburst and 
thus drain even more energy from the shock-heated mat- 
ter USEil. 

While v losses damp the shock in this early phase, the 
situation changes some 50 ms later. As more stellar mat- 
ter falls onto the collapsed inner core, the shock is pushed 
to larger radii and the density and temperature behind 
the shock decrease. On the other hand, the central core 
begins to settle and heats up, thus radiating more ener- 
getic neutrinos. Both effects lead to the situation that 
v e and v e can now be absorbed with a small probability 
(10-20%) by free neutrons and protons behind the shock. 
A region of v heating between the so-called "gain radius" 
and the shock front develops Q ■ If the v energy deposi- 
tion is efficient enough, the stalled shock can be revived 
and drives a "delayed" explosion. 

The success of pioneering calculations could be 



maintained in later models only by invoking neutron- 
finger convective instabilities in the PNS |3J|]. These 
boost the v luminosities and thus enhance the z/-energy 
transfer to the shock. Explosion energies similar to those 
of observed SNe required even higher v fluxes. It was pro- 
posed that these could be obtained when pions appear in 
large numbers in the PNS matter The existence 

of neutron-finger instabilities, however, depends on very 
specific thermodynamical properties of the equation of 
state and on the details of the v transport The for- 
mation of pions in hot PNS matter, on the other hand, 
is highly uncertain and requires particular assumptions 
about their dispersion relation. 

While all simulations addressed so far were performed 
in one dimension (ID) assuming spherical symmetry 
(neutron-finger convection was treated by a mixing- 
length approach), SN 1987A provided evidence for large- 
scale mixing processes which carried radioactive nuclei 
from the region of their formation near the PNS into 
the helium and hydrogen shells of the exploding star. 
Simulations suggested that their origin may be linked 
to hydrodynamic instabilities behind the stalled shock 
already during the first second of the explosion [lfij . 
Two-dimensional (2D) || and most recently also three- 
dimensional (3D) ^3 models that take into account v 
effects then showed that convective overturn indeed de- 
velops in the v heating region and is helpful for shock 
revival, thus making explosions possible even when spher- 
ically symmetric models fail 19j. 

The multi-dimensional situation is generically differ- 
ent because it allows accretion to continue while shock 
expansion already sets in. Narrow downflows bring cold, 
low-entropy matter close to the gain radius, where the v 
energy deposition is strongest. At the same time heated 
matter can rise in buoyant bubbles, thus pushing the 
shock farther out and reducing energy loss by the ree- 
mission of neutrinos. Although this increases the effi- 
ciency of ^-energy transfer, convection is still no guar- 
antee that explosions occur [19j . A particular concern 
with all multi-dimensional models which yielded explo- 
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sions was the much simplified treatment of the v trans- 
port by grey diffusion schemes, which was inferior to the 
more elaborate multi-group tr ansp ort description used in 
unsuccessful spherical models |28] . 

Recently it has become possible to go a step further 
and solve the time-dependent Boltzmann equation for v 
transport in ID hydrodynamic simulations with Newto- 
nian [22I HH EH and relativistic gravity ,24] . It turned 
out that even this improvement does not lead to explo- 
sions. The question, however, remained whether con- 
vection might bring the models to a success. Here we 
present the first 2D simulations which were performed 
with a Boltzmann solver for the v transport. At the same 
time we have also upgraded the description of ^-matter 
interactions compared to the conventional treatment of 

Refs. HHim. 

Numerical techniques and input physics. For the in- 
tegration of the equations of hydrodynamics we employ 
the Newtonian finite-volume code PROMETHEUS [Hi 
This second-order, time-explicit Godunov scheme is a di- 
rect Eulerian implementation of the Piecewise Parabolic 
Method (PPM) and is based on a Riemann solver. 

The Boltzmann solver scheme for v and v of all 
three flavors is described in detail in Ref. [3l| . In 
multi-dimensional simulations in spherical coordinates, 
we solve for each latitude 9 of the numerical grid the 
monochromatic moment equations for the radial trans- 
port of v number, energy, and momentum. This set of 
equations is closed by a variable Eddington factor that 
is calculated from the solution of the Boltzmann equa- 
tion on an angularly averaged stellar background. In 
addition, we had to go an important step beyond this 
"ray-by-ray" approximation of multi-dimensional trans- 
port. Physical constraints, namely the conservation of 
lcpton number and entropy within adiabatically moving 
fluid elements and the stability of regions which should 
not develop convection according to a stability analysis, 
made it necessary to take into account the coupling of 
neighboring rays at least by lateral advection terms and 
v pressure gradients Q. 

General relativistic effects are treated approximately 
by modifying the gravitational potential with correction 
terms due to pressure and energy of the stellar medium 
and neutrinos, which are deduced from a comparison of 
the Newtonian and relativistic equations of motion |3lj |. 
The v transport contains gravitational redshift and time 
dilation, but ignores the distinction between coordinate 
radius and proper radius. This is necessary for coupling 
the transport code to our basically Newtonian hydrody- 
namics. Comparison with fully relativistic ID simula- 
tions showed that these approximations work well at least 
when the deviations of the metric coefficients from unity 
are moderate [23). 

Improving the description of v interactions Q com- 
pared to the "standard" opacities 0, |2(| we added 
j/-pair creation (and annihilation) by nucleon-nucleon 



bremsstrahlung |14j . scattering of v^, v^, v T , and v T off 
v e and z> e , and pair annihilation reactions between neutri- 
nos of different flavors . We also take into account the 
detailed reaction kinematics with nucleon thermal mo- 
tions, recoil, and fermion phase-space blocking effects in 
the charged- and neutral-current ^-nucleon interactions. 
Moreover, nuclear correlations ||, weak-magnetism cor- 
rections |l7j, an d the reduction of the nucleon effective 
mass and possible quenching of the axial- vector coupling 
in nuclear matter [lflj are included, too. 

TABLE I: Parameters of computed 2D models for progenitor 
stars with different masses. Mf c is the iron core mass, Msi+o 
the mass interior to the inner edge of the oxygen-rich Si-shell, 
Hi the angular velocity of the Fe-core prior to collapse, 80 and 
9i are the polar angles of the lateral grid boundaries, and Ne 
is the number of grid points in lateral direction. 
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15 


1.28 


1.43 
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64 
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20 


1.46 


1.46 
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Our current SN models are calculated with the equa- 
tion of state (EoS) of Ref. using an incompressibility 
of bulk nuclear matter of 180 MeV (other values make 
only minor differences 32]). This EoS is based on a 
compressible liquid-drop model including nuclei, nucle- 
ons, e~ , e + , and photons. 

The 2D simulations were performed with a spherical 
coordinate grid with about 400 radial zones and assuming 
azimuthal symmetry around the polar axis. A wedge of 
roughly 90 degrees was chosen in the angular direction 
between the angles 60 and 0\ and with the number of 
mesh points, Ng, as given in Tabled For the nonrotating 
models periodic boundary conditions were taken. For the 
rotating model reflecting boundaries at the axis and the 
equator were imposed. 

We considered three (solar metallicity) progenitors 
with main sequence masses of 11.2 M©, 15 M©, and 
20 M© [34| . The model specific parameters are given in 
Table [I] For the transport we used an energy grid of 
17 geometrically spaced bins with centers from 2 MeV 
to 333 MeV. Higher energy resolution was tested without 
giving significant differences. Model sl5p is a calcula- 
tion where random density perturbations with 1% am- 
plitude were imposed on the pre-collapse core to follow 
their growth during infall |l8|. In the rotating Model 
sl5r the angular frequency was assumed to be constant 
in the Fe and Si core and decreasing like r -3 / 2 outside of 
1750km (1.43 M©). These choices and the adopted ro- 
tation rate (cf. Table [I| are basically in agreement with 
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predictions from stellar evolution models [l 
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FIG. 1: Snapshots of the stellar structure for the rotating 
Model sl5r at 198 ms after shock formation. The left panels 
show the rotational velocity (top) and the fluctuations of en- 
tropy (in per cent) versus the enclosed mass, emphasizing the 
conditions inside the neutron star. The right panels display 
the rotational velocity (top) and the entropy (in Icb per nu- 
cleon) as functions of radius. The arrows indicate the velocity 
field, the white line marks the shock front. 
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FIG. 2: Shock radii R s (averaged over polar angles) vs. post- 
bounce time i p b. The 2D models (bold lines) are compared 
to the corresponding ID simulations (thin lines). 

Results. Convective overturn starts to appear in the 
region of v heating behind the shock about 25 ms after 
shock formation, and begins to affect the SN dynamics 
about 50 ms later. Expanding bubbles of heated mat- 
ter deform the shock (Fig. ^) and lead to a transient 
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FIG. 3: Luminosities and mean energies (defined as ratio of 
energy flux to number flux) for v e , i/ e and ^.t, v^, t vs. time 
for Models sll.2 (top) and sl5r (measured by an observer 
comoving with the stellar medium at 500 km) . The left panels 
(left scale) show the prompt v £ burst, the right panels enlarge 
the post-bounce evolution. The thin lines represent results of 
ID simulations for comparison. 



increase of the average shock radius (Fig. |2J). The dif- 
ference to spherically symmetric simulations, however, is 
small in case of Models sl5 and s20 because the con- 
vectively unstable layer between gain radius and shock 
remains narrow and the overturn motions never become 
very strong. The largest relative change occurs when the 
big entropy discontinuity at the inner boundary of the 
oxygen enriched shell (Msi+o; Table [IJ of the progenitor 
star crosses the shock. The sudden decrease of the den- 
sity and ram pressure of the infalling matter allows the 
shock to transiently expand also in ID simulations. In 
Model sll.2 the corresponding growth of the shock ra- 
dius is nearly 50 km. This in turn widens the gain layer 
and thus allows convection to strengthen. The combined 
effects lead to a much larger shock radius than in the ID 
case (Fig. [2J • A similar result is caused by centrifugal 
forces in the rotating Model sl5r (Fig. [J). Here decel- 
erated infall of matter near the equatorial plane and a 
reduced density near the poles create a more favorable 
situation for violent convection (Fig. ^| . 

Ledoux convection sets in below the neutrinosphere al- 
ready about 40 ms after bounce. It is persistent until the 
end of our simulations and slowly digs farther into the 
star (cf. Fig. ^ lower left panel) in agreement with pre- 
vious 2D simulations [20j . But the active layer is rather 
deep inside the PNS (at a density above 10 12 gem -3 ) and 
is surrounded by a convectively stable shell in which the 
surface fluxes of v e and v e are mainly built up. Therefore 
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the PNS convection has little influence on the emission 
of these neutrinos and is irrelevant for the SN dynam- 
ics. The differences in the v e and D e luminosities and 
mean energies relative to the ID case (Fig.[3J) are caused 
by rotation and matter downflows from the shock to the 
neutrinospheric region. 

Random density perturbations in the progenitor 
(Model sl5p) increase in the supersonically infalling lay- 
ers and are damped in the subsonically collapsing inner 
core as predicted by linear analysis 01 • We have not dis- 
covered any influence on the development of convection. 
Also the shock radii of Models sl5p and sl5 are nearly 
identical. Higher angular resolution allows convection to 
grow somewhat faster because of reduced numerical vis- 
cosity. Within the limits of our tests the later evolution 
in cases of strong convection exhibits quantitative but no 
qualitative differences (see Models sll.2 in Fig.^J). 

Conclusions. Our simulations show that an accurate 
treatment of the neutrino physics does not yield suffi- 
ciently efficient v-eneigy transfer behind the stalled SN 
shock to produce explosions, even though convection oc- 
curs below the neutrinosphere and in the ^-heating re- 
gion. This finding confirms concerns 28] that the success 
of previous models was favored by gross simplifications 
of the v transport. 

None of the included effects is therefore able to cause 
explosions. Truly multi-dimensional transport Q or full 
relativity [13 are not likely to change the situation. But 
the long-time evolution of the shock is sensitive to the v 
emission that originates from the neutrinospheric layer. 
The EoS properties and v interactions at densities below 
some 10 13 gem -3 are therefore particularly relevant. The 
structure, temperature, and convective stability of this 
layer depend also on the compactness of the PNS and 
thus on the uncertain physics in its supranuclear core. 
The influence of the EoS on the post-bounce dynamics, 
however, has not been explored extensively so far. Neces- 
sary improvements of weak interactions on nuclei include 
electron capture rates [2lJ and inelastic neutrino scat- 
tering. Moreover, generically 3D phenomena cannot be 
studied with 2D hydrodynamics. Also the properties of 
the progenitor cores need reexamination. 

The current paradigm for explaining massive star ex- 
plosions would have to be revised if the ^-driven mecha- 
nism were fundamentally flawed. SNe might be aided by, 
e.g., magnetohydrodynamic processes jl| or even more 
exotic physics. Although it is still too early for this con- 
clusion, such investigations deserve more interest. 

We are indebted to K. Takahashi and C. Horowitz 
for help in improving ^-nucleon interactions. Support 
by the Sondcrforschungsbereich 375 on "Astroparticle 
Physics" of the Deutsche Forschungsgemeinschaft is ac- 
knowledged. The simulations were done on the IBM "Re- 
gatta" supercomputer of the Rechenzentrum Garching. 
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